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COMPUTATION OF SUB-SATELLITE POINTS 
FROM ORBITAL ELEMENTS 
By Richai'd H. Chi'ist 
John F. Kennedy Space Center 

SUMMARY 

This technical note presents the computer program used by the Computation Branch at 
the John F. Kennedy Space Center (KSC), NASA for computing predicted sub-satellite 
points of radiating satellites for launch interference purposes. 

INTRODUCTION 

Frequently, it is necessary to predict launch interference from radiating satellites. 

The predicted launch interference information is desired in the form of sub-satellite points 
and look angles. This technical note describes the method used by the Computation Branch 
of KSC for computing the sub-satellite points from Prediction Space Elements provided by 
the Goddard Space Flight Center (GSFC), NASA. 

The orbital elements at the time of interest, t, are interpolated from the given ephe- 
meris; the orbital elements at t are transformed to a position vector at t; and the sub-satel- 
lite points at t are computed. The output of the computer program may be then used for 
standard look angle computations. 

Thanks are due to Mr. T. P. Gorman, Chief of the Advanced Orbital Programming 
Branch, GSFC for providing information on the prediction elements; and to Mr. W. N. 
Weston of GSFC who provided check data. 

SYMBOLS 


a semimajor axis of orbit 

e eccentricity 

i inclination 

Q right ascension of ascending node 

tu argument of perigee 

M mean anomaly 

E eccentric anomaly 

P period (anomalistic) 

P period derivative 

r^ radius vector 

x,y,z geocentric, inertial, right-handed, orthogonal position components; z is co- 
incident with the polar axis and x is directed toward the vernal equinox 
J_,_j,l^ unit vectors lying along the x,y,z coordinate axes, respectively 

e eccentricity of reference spheroid 



$' geocentric latitude 

$ geodetic latitude 

X longitude 

h height above spheroid 

^ rotational rate of earth 

EL general term referring to orbital elements 

t time of interest 

tj ,tg times of epochs 1, 2, and 3 respectively 

R .A. right ascension 

S.T. sidereal time 

tp reference time 

U.T. universal time 

C.U.L. canonical unit of length = 6378. 165 km 
J . D . S . J ul ian Date for S pace 

J.D.S. k J.D. - 2,436,099.5 days 
All times are expressed in terms of J.D.S. 
unless specified otherwise 


COMPUTER PROGRAM EQUATIONS 
The following are programmed digital computer equations; 

Subprograms ELIN and EVERET 
Interpolation for aft), eft), i(t),Q ft), andouft) - 



ft) 

= a(t). 

Ekx 

= ^1 

^ Ek, 

= ^3 

/Eka 

LU 

ft) 

= e(t). 

Ekx 

= ^1 

/Ek, 

= 63 

/Eka 

EL 3 

ft) 

= i(t). 

Ekx 

“ *1 

Ek. 


'EL33 

Ek 

(t) 

=^(t). 

Eki 

= k 

, EL^3 


/ EL 43 

Ek 

(t) 

=«J(t), 

Ekx 


,EL33 


/ ^1-53 


R = 1 - S 


2 



A® =EL. - 2EL. +EL. 

0 IS IS 1 1- 


EL| (0 = R + S EL, , 

S (S= - 1) S 

A 

+ 3 ! 0 


where i = 1, 2, • • v 6. 


Computation of M(t) - 


Subprogram AMIN 

Cl = 518,400 
Cs =373.248 


For the interval between t. and t. 

I I 4 1 


C P- 

M (t) = Mj 4-^ (t - tj) - Cs (t - tj)® 

rj Kj 


Subprogram EKEP 

Solution of Kepler's equation - 

e (t) sin M (t) 

y e2 (t) 4 1 - 2e (t) cos M (t) 

E (t), = M (t) 4 Z - ^ cot M (t) 

^ 6 

M (t) 4 e (t) sin E (t). - E (t). 

E (t). , = E (t). 4 ! •- 

• + ^ ' 1 - e (t) cos E (t) j 
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Continue iterating until |E(t)j_j.^ -E(t)jj < 0 .2 x 10~’ rad. 

Subprogram PVOE 

Components of position vector - 

X (t) = a (t) I [cos E (t) - e (t)] [cos c«(t) cosO (t) 

- sin ou(t) Sinn (t) cos i (t) ] 

+ [ 1 - e® (t)]^ sin E (t) [- sinou (t) cosO (t) 

- cosuj (t) sin n(t) cos i (t) ]| 

y (t) = a (t) I [cos E (t) - e (t)] [cos («(t) sinO (t) 

+ sincu(t) cos n(t) cos i (t) ] 

+ [ 1 - e® (t)]^ sin E (t) [ - sinou (t) sin Q(t) 

+ cos uu(t) cos n(t) cos i (t) ]j. 

z (t) = a (t) |[cos E (t) - e (t)] [sinoo (t) sin i (t) ] 

+ [1 - e® (t)]^ sin E (t) [cos uu (t) sin i (t) ] j. 

Subprogram GEODT 

Sub-satellite points - 

_i y (t) 

R . A . (t) = tan 

x(t) 

Atp = t - tp 

Xq = S. T. atO^ U. T. 

X(t) = R. A. (t) - X^ 

r (t) =yx®(t) + y®(t) +z®(t) 
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I 





where 




— A , - A 
0 +5 “2 

Everett's Second Central Difference formula is given as 


V. - DY a- R(R^ - ) A® , R(R^ - f ) (R^ - f ) a" 

Xj ^Xq I 0 1 0 


3! 


5! 


+ SX^ + 


S(SP - f ) A 


3! 


, S(S - 1 ) (S - 2 ) 

1 + 

5! 


where S = — ' ° , R = 1 - S , 

W 


and W is the time interval between entries in the table. 


( 1 ) 


For our application, only three X's are available; therefore, the terms containing the 
fourth differences may be eliminated. Letting a®q = ,and using the notation intro- 
duced in the previous section (Computer Program Equations), we have 


EL (t) = R ELg + 


R(R® - 1) 
3! 


,+ S EL, + 


S(S® -1) 


3! 


^0 < 2 ) 


where 


^ Ao-i - i-A 


0 - +5 


= X, -2Xo + X_, 

A® =ELs - 2 ELs +ELi 
0 


This scheme was compared with a second degree least squares fit and agreement to 
0.5 X 10*’^ degree was reached for Q(t). 


Computation of M (t) 

The mean anomaly at t, M (t), assuming a constant period, P, is determined from 
Kepler's equation 
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(3) 


M (t) = Mj (t - tj) 


Since P is available as input, equation (3) may be expanded by Taylor's series 


f (x) — f (a) + f (a) (x - a) H — ^ f " (a) (x - a)' 
+ • • f (a) (x - a) " 


(4) 


f (x) = M (t) 


2n 


f (a) = Mj d — pT“(t - tj) 

p (t - tj) 


= Mj 


(5) 


If (t - t|) is in days, and Pj in minutes; 

= 57.2957795 x 1440 x 2 n = 518,400 


f (a) 

f (a) (x - a) 


Pj Pj^ 

= |^(t - tif 

C, 2 C« • 

^1 '^\ 


( 6 ) 


If (t - t|) Is In days , Pj In minutes, Pj in mlcrodays per day; 


C =TTx 57.2957795 x (1440f X 10’® 
= 373.248 
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(7) 


Substitution of equations (5) and (6) into equation (4) yields the desired expression 

P* 

M (t) = Mj + — - tj) - Cs (t - tj)® 

Equation ( 7 ) is the expression given in reference 1. 


Solution of Kepler's Equation 

The number of methods devised for the solution of Kepler's equation exceeds one 
hundred. The method presented below is an iterative scheme well suited fora high- 
speed digital computer. 


M (t) = E (t) - e (t) sin E (t) 
It is desired to solve equation ( 8 ) for E (t). 

The Newton-Raphson formula is given as 

f (Xj) 

' f (Xi) 


(8) 


(9) 


Application of equation (9) to equation (8) gives 

M (t) + e (t) sin E (t)j - E (t),- 


=E(t)i + 


1 - e (t) cos E (t) 


Encke's first approximation is very good, seen as 


(10) 


E (t), = M (t) Z - T cot M (t) 
^ 6 

e (t) sin M (t) 


where 


Z = 


yes (t) + 1 - 2 e (t) cos M (t) 


Continue iterating equation (10) until 
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<0.2 X 10“^ rad. 


Position Vector at t 

The following equations (11 - 15) are taken from reference 1 . 

Q (t) = COS n(t) + sin n(t) J_ 
g (t) = cos i (t) + [ sin i (t)] [_Q^ (t) x k ] 

_g_(t) = cosuu (t) Q(t) + [sinuu(t)]|g (t) xQ(t)] 
jg_ (t) = _^ (t) X _g_ (t) 

_r_ (t) = a (t) I [cos E (t) - e (t)] (t) 


- es (t) sin E (t)_g_(t)} 


( 11 ) 

(12) 

(13) 

(14) 

(15) 


Performing the indicated operations in equations (11-14) and substituting in equation 
(15) gives the desired expression for_r_ (t) as 

_r_ (t) = a (t) I [cos E (t) - e (t) ] [cosiw (t) cos Q (t) 

- sin (ju (t) sin Q (t) cos i (t) ] 

+ [1 - e® (t)]^ sin E (t) [ - sinuu (t) cos n (t) 

- cosuu (t) sino(t) cos i (t)]| 

+ j a (t) |[cos E (t) - e (t) ] [cosuu (t) sinn (t) 

+ sinuu (t) cos Q (t) cos i (t) ] 

+ [1 - e®(t)]^ sin E (t) [- sinuu(t) sin n(t) 

+ cos uu (t) cosn (t) cos i (t)] I 
+ k a (t) |[cos E (t) - e (t) ] [sin uu (t) sin i (t) ] 


+[1 (t)]^sinE (t) [ cos uj (Osin i(t)]}(16) 
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Sub-satellite Points 


The method used to compute geocentric latitude, and longitude, X, is similar to 
that described in reference 2 , as 


$ ' (t) = tan 


- X 


2 (t) 


7x2 (t) + y2 (t) 


(17) 


R. A. (t) = tan“' 

X (t) 


(18) 


Xq = S. T. at 0" U, T. 


(19) 


then 

where 


X(t) = R. A. (t) - Xq -(J) Atp 


AtR=t-tR 


( 20 ) 


For the determination of geodetic latitude,#, the first two terms of the non-iterative 
scheme presented in reference 3 are used as 


where 


§(t) = #'(t) + ag (t) sin 2#'(t) + a 4(0 sin 4#'(t) 

+ 128 e^ + 60 e® + 35 


g ^ 512 ^^ •a c As i 


1024r (t) 

^ fg® A8 ^ 


32r2 (t) 
-1 




“ 


1024r(t) 


256r3 (t) 

[64 + 48 +35 e® ] 


4 ^®+ 3 ^"], 


+ [4 ^^ + 2 +e® j 


16r2 (t) 

+ 


15 ^® 


gs 


256rMt) 16r(t) 


( 21 ) 
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and 


r (t) + y^(t) + z^(t) , in C. U. L. 

Finally for height above spheroid, h 

sin $ (t) 1 - e^ 

h (t) = r (t) ^ 

sin $ (t) [1 -^^sin^ $(t )]2 


( 22 ) 


COMPUTER PROGRAM OPERATING INSTRUCTIONS 

Program Description 

Program Identification - SPOE 1 
Computer - GE 235 
Program Library - 111141 
Type of Coding - FORTRAN II 


Program Input (From Cards) 


Card No. 

Columns 

Information 

Mode 

Remarks 

1 

34 

- 45 

Date 

4A3 

Information on cards 






1 and 2 for identifi 

2 

27 

- 50 

Satellite ID 

8A3 

cation only 

3 

1 

- 3 

tR (y.) 

13 

Reference Time 


4 

- 6 

tj^ (mo.) 

13 



7 

- 9 

tR (dy.) 

13 



10 

- 12 

tp (hr.) 

13 



13 

- 15 

tp (min.) 

13 



16 

- 22 

tp (sec.) 

F7.3 



23 

- 25 

S. T. (hr.) 

13 

S. T. at O'^U. T. 


26 

- 28 

S . T. (min.) 

13 

(Reference U . T .) 


29 

- 35 

S . T. (sec.) 

F7.3 



36 

- 46 

*^R 

F11.5 

Reference time in 






J. D. S. 


47 

- 61 

At 

E15.9 

Time increment 






desired in days 
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Program Input (From Cards) (Cont'd) 


Card No. Columns 

4 1-15 


16 - 30 


5 1-11 
12 - 22 
23 - 33 

34 - 44 
45 - 55 
56 - 66 

6 1-11 
12 - 22 
23 - 33 

7 1-11 
12 - 22 
23 - 33 

34 - 44 
45 - 55 
56 - 66 

8 1-11 
12 - 22 
23 - 33 

34 - 44 
45 - 55 
56 - 66 

9 1-11 

12 - 22 
23 - 33 


Information 

Mode 

A 

UU 

E15.9 

r 

E15.9 

ti 

F11.5 


F11.5 

ts 

F11.5 

Pi 

F11.5 

K 

F11.5 

Ps 

F11.5 

p. 

F11.5 

Ps 

F11.5 

P3 

F11.5 


F11.8 

^2 

F11.8 

as 

F11.8 


F11.8 


F11.8 

^3 

F11.8 


F11.6 

*2 

F11.6 

is 

F11.6 


F11.6 

fig 

F11.6 

Q3 

F11.6 

U)l 

F11.6 

0)2 

F11.6 

OJ3 

F11.6 


Remarks 

Rotational rate of earth 
in radians per second 

z 

(Eccentricity) of 
spheroid 

Epoch times in J. D. S. 


Period in minutes 


Period derivative in 
microdays per day 

Semimajor axis of orbit 
in earth radii 


Eccentricity of orbit 


Inclination in degrees 


Right ascension of 
ascending node in 
degrees 

Argument of perigee 
in degrees 
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Program Input (From Cards) (Cont'd) 


Card No. 

Columns 

Information 

Mode 

Remarks 


34 - 44 

Ml 

F11.6 

Mean anomaly in 


45 - 55 


F11.6 

degrees 


56 - 66 

M3 

F11.6 


10 

1-3 

CNTRL 

F3.0 

Output control card 
-1. print only 


0. write tape only 
+ 1. print and write 
tape 


Program Output 

Option 1: Print only. 

Option 2: Write tape only; output tape on plug 2, unit 1. 

Option 3: Print and write tape; output tape on plug 2, unit 1. 

Card 10 specifies output option. Output record is t, § (t), X (t), and h (t). Input data 
are printed on all options. 

Sample Test Case 

Input. - The following listing is a computer printout of a sample input giving the 
Prediction Space Ejements supplied by GSFC for program input; 

DEC 28 1964 

SATELLITE RELAY 2 


64 12 01 00 00 00,000 04 39 30.641 2631.00000 . 6OOOQQO0QE-Q2 


.729211590E 

-04 .669342162E-02 




2631.00000 

2639.00000 

2647.00000 

0194.71578 

0194.71538 

0194.71497 

0000.03520- 

0000.03537- 

0000.03559 




1.74488160 

1.74487910 

1.74487650 

0.23957545 

0.23950386 

0.23943548 

046.326254 

046.327234 

046.32816^ 

236.668682 

227,838464 

219.008515 

172.065590 

180.909426 

189.755772 

318.158721 

016.916503 

075.718748 


+ 1 
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Output. - The following listing is a sample computer printout of sub-satellite points. 


JOHN F. KENNEDY SPACE CENTER 

computation branch 

SUB-SATELLfTE POINTS 
DEC 2a 1964 

SATELLITE RELAY 2 


REFERENCE TIME 

YR 

CAL DAT UT2W 64 

MO 

12 

DY HR MM SS.SSS 
1 0 0 0 . 

S.T. at 0 hr U.T. 
HR MM SS.SSS 
4 39 30.641 


J.D.S. UT2W 2631. 

000 00 



PREDICTION SPACE ELEMENTS FROM GODDARD 



EPOCH 



T-ONE 

T-TWO 

t-three 

J.D.S. 

UT2W 


2631.00000 

2639.00000 

2647.00000 

PERIOD 

MIN 


194.71578 

194.71538 

194 , 71497 

PERIOD DER 

MD/D 


-0.03520 

-0.03537 

-0.03559 

ECCENTRICITY 



0.23957545 

0.23950386 

0.23943548 

INCLINATION 

DEG 


46.326254 

46.327234 „ 

46.328169 

RA ASC NODE 

DEG 


236.668682 

227.838464 

219.008515 

ARG PERIGEE 

DEG 


172.065590 

leO .909426 

189.755772 

MEAN ANOMALY 

DEG 


318.158721 

16.916503 

75.718748 

SEMIMAJ AXIS 

ER 


1.74488160 

1.74487910 

.1 .74487650 
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TIME[ J.D.S. 3 

LATCDEG] 

LON[nEQ] 

HEIGHT [KM] 


2631 . OOUUU 

43.78 

- 78.84 

3148.1 


2631 .00600 

33,71 

- 54.73 

2532.7 


2631 . 0119 V 

17.52 

-^.0 0 

2155.0 


2631 . 0179 V 

- 1.62 

- 18.18 

2110.4 


2631 . 02399 

- 20.03 

— 1.63 

2412.3 


2631 .02998 

- 34.64 

16.99 

2981.5 


2631 . 03598 

- 43.59 

38.58 

3699.8 


2631 ,04198 

- 46.44 

61.06 

4462.9 


2631.04797 

- 44.44 

80.56 

5198.0 


2631 . 0539 / 

- 39.56 

95.45 

5859.1 


2631 . 0599 / 

- 33.30 

106.40 

6418.9 


2631 . 06596 

- 26.45 

114.65 

6861.8 


2631.07196 

- 19.36 

121.20 

7178.5 


2631.07796 

- 12.18 

126.71 

7364.3 


2631 .08395 

- 4.92 

131.70 

7416.7 


2631 . 08995 

2.43 

136.54 

7335.2 


2631 . 09595 

9.91 

141.62 

7120.8 


2631.10194 

17.52 

147.36 

6776.5 


2631.10794 

25.24 

154.36 

6307.9 


2631.11394 

32.90 

l 63 .53 

5725.4 


2631.11995 

39.97 

176.25 

5047 . 0 


2631.12593 

45.18 

- 165.69 

4303.0 


2631.13193 

46.08 

- 141 .85 

3544.0 


2631 . 13792 

40.08 

- 116.26 

2848.7 


2631.14392 

26.91 

- 9 4 . 1 . 5 

2325.9 


2631.14992 

8.92 

- 76.07 

2090.4 

— 
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TIME[ J.D.S. J 

LaT [DEG] 

LGN[DER] 

HElGHTtKM] 

2631.15591 

-10 .33 

-59.69 

2206.7 

2631.16191 

-27.33 

-42.43 

2642.9 

2631.16791 

-39.51 

-22.39 

3294.0 

2631.1739U 

-45.62 

0. j 0 

4043.6 

2631.1799U 

-46.00 

21.65 

4801.2 

2631.18590 

-42.43 

39.07 

5507.3 

2631.19189 

-36.77 

51.99 

6125.2 

2631 . 19789 

-30.16 

61.56 

6633*8 

2631 .20389 

-23.17”^ 

68.93 

7 0 2 0.6 

2631.20988 

-16.03 

74 . 93 

7278.7 

2631.21588 

-8.80 

80.15 

7404.5 

2631.22188 

-1.51 

85 .02 

7396.5 

2631.2278/ 

5.89 

89.92 

7254.8 

2631.2338/ 

13.43 

95.24 

6981.3 

2631.2398/ 

21.10 

101.48 

6580 . 0 

2631 . 24580 

28.82 

109.33 

6058.4 

2631.25180 

36.30 

119.92 

5429.8 

2631.25780 

4 2.72 

134.87 

4717.0 

2631.26385 

46.31 

155.73 

3950.5 

2631,26980 

4 4.28 

-178.90 

3216.6 

2631.27580 

34.85 

- 154.51 

2584.3 

2631.28185 

19.09 

- 134.45 

2178.3 

2631.28784 

0.06 

-117.49 

2098.9 

2631.29384 

^18756 

-101 . 00 

2369.2 

2631.29984 

-33.61 

-82.60 

2917.3 

2631.30583 

-43.09 

-61.17 

3626.1 
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TIMEU.D.S. j 

LaTLDEG] 

LON [DEG] 

HEIGHT[KM) 

26 ^ 1 * ■ 3 Xl . dN 5 

- 46.41 

- 38.56 

4388 . 6 

2631.31763 

- 44.72 

- 18 .66 

5128.7 

2631.32382 

- 39.98 

- 3.39 

5798.4 

2631 . 3298 ^ 

- 33.78 

7.84 

6369 . 0 

2631 .33582 

- 26.94 

16.27 

6823.8 

2631.34181 

- 19.86 

22.92 

7153.2 

2631.34781 

- 12.67 

28.50 

7351.9 

2631 .35381 

- 5.41 

33.51 

7417.4 

2631.35980 

1.93 

36.34 

7349 . 1 

2631 .36580 

9.38 

43.39 

7147.6 

2631.37180 

16.98 

49 . 06 

6815.8 

2631.37779 

24.68 

55.93 

6359 . 0 

2631 .38379 

32.34 

64.87 

5787.1 

2631.38979 

39,47 

77.22 

5116.9 

2631.39578 

44.88 

94.76 

4377.5 

2631,40178 

46.24 

118.15 

3616.9 

2631 .40778 

40.89 

143.77 

2911.1 

2631.41377 

28.27 

166.24 

2366.0 

2631.41977 

10.58 

- 175.43 

2098.3 

2631.42577 

- 8.71 

- 159 .01 

2180 . 0 

2631.43178 

- 26.04 

- 141 .89 

2588.8 

2631 . 43776 

- 38.72 

- 122.07 

3224.2 

2631 . 4437 b 

- 45.35 

- 99.63 

3968.7 

2631.44975 

- 46.14 

- 77.80 

4728.6 

2631.45575 

- 42.80 

- 59.97 

5441.5 

2631.46175 

- 37.23 

- 46.70 

6069.3 
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COMPUTER PROGRAM FLOW CHART 


Main Program SPOEl 



• COMPUTE I 
a(t), e(t), 

i(t), m. & 

o>(t) 

USING 

ELIN & EVERET 


1. COMPUTE 
M(t) USING 

AMIN 

2. CONVERT 
M(t) TO 
RADIANS 


COMPUTE 
E(t) USING 

EKEP 


CONVERT 
i(t), fZ(t) and o>(t) 
TO RADIANS 


COMPUTE 

x(t), y(t), z(t) 

USING PVOE 


COMPUTE 
^{♦)/ M*) °nd h(t) 
USING GEODT 


1. CONVERT 0{t) 
ond A(t) to deg. 

2. CONVERT h(t) to Icm 


test 


m 

CNTRL 



+1 

V 
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3 


WRITE BINARY 
TAPE 

t, A(t), h(t) 




Subprogram ELIN 


o(t)= ELJt) 
e(t)=EL^(t) 
i(t)=EL^(t) 
n(t)= EL^(t) 
a>(t)= ELJt) 


Subprogram EVERET 










Subprogram AMIN 





P_ = P , P^ = P^ 



P, = P^, P^ = P^ 



COMPUTE 

M(t) 




RETURN TO 
\MAIN PROGRANy 





Subprogram EKEP 



2] 
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COMPUTER PROGRAM LISTING 


* 

c 

c 

c 

_c 

c 

c 

c 

c 

c 


c 


FORTRAN 

SPOEl- COMPUTATION OF SUB-SATELLITE POINTS FROM ORBITAL ELEMENTS 

PROGRAMMED BY RICHARD H» CHRIST 

UTILIZES - FUNCTION AMIN 
FUNCTION EKEP 
SUBROUTINE ELIN 

FUNCTION EVERET 

SUBROUTINE PVOE 

SUBROUTINE GEODT 

DIMENSION EL(6»3) *AL(5»3) *A1(4) »A2(8) 

READ INPUT DATA 

READ 123*A1 

READ 124»A2 

READ 100» IYR» IMO» IDY* IHR. IMM»SS» I SHR ♦ I SMM . SSS ♦ RT . DELTAT 

READ 101»VE.ES 

READ 102.T1.T2.T3.P1.P2.P3 

READ 102.D1»D2.D3 

READ 103. ( ( EL( I . J) . J=1.3 ) *1 = 1*2) 

READ 104*( {EL( I *J) *J = 1*3) *1=3*6) 

READ 129.CNTRL 

PRINT INPUT DATA 

PRINT 105 

PRINT 121 

PRINT 122 

PRINT 127 

PRINT 123.A1 

PRINT 124.A2 _ 

PRINT 125 
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I 


PRINT 106 


PRINT 107 

PRINT 108.IYR.IMO»IDY.IHR.IMM.SS . I SHR . I SMM » SSS 

PRINT 120. RT 

PRINT 109 

PRINT 110 

PRINT 111.T1.T2.T3 

PRINT 112.P1.P2.P3 

PRINT 113.D1.D2.D3 

PRINT 114.(EL(2»J) .J=1.3) 

PRINT 115.(EL(3.J) .J = 1.3) 

^INT 116.(EL(4.J) .J = 1.3) 

PRINT 117.(EL(5.J) .J=l»3) 

^PRINT 118. ( EL ( 6 . J ) . J = 1.3 ) 

PRINT 119.(EL( l.J) .J = 1.3) 

C ^FORMAT STATEMENTS TO READ AND PRINT INPUT 

100 FORMAT (3I3»2(2I3.F7.3) .F11.5.E15.9) 

101 FORMAT (2E15.9) 

1^0 2 FORMAT (6F11.5) 

^3 FORMAT (6F11.8) 

lOA FORMAT (6F11.6) 

1^5 FORMAT (IHI) 

106 FORMAT ( 58HREFERENCE TIME S.T. AT 0 HR U 

l.T. ) 

107 F ORMAT (16X.38H YR MO DY HR MM SS.SSS HR MM SS.SSS) 

^8 FORMAT (16H CAL DAT UT2W . 5 I 3 , F7. 3 . 3X . 2 1 3 . F7 . 3 ) 

109 FORMAT ( / / / 38HPRED I CT I ON SPACE ELEMENTS FROM GODDARD) 

110 FORMAT (/72HEPOCH T-ONE T-TWO 
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1 


T-THREE) 


111 

FORMAT 

(/17HJ.D.S. 

UT2W.4X*3(F11.5*10X 

) ) 

112 

FORMAT 

( /17HPERIOD 

MIN*4X*3(F11.5*10X 

) ) 

113 

FORMAT 

(/17HPERI00 DER 

MD/D.4X*3(F11.'5*10X 

) ) 

11^ 

FORMAT 

(/17HECCENTRICITY 

.7X.3(F11.8»10X 

) ) 

115 

FORMAT 

( /17HINCLINATION 

DEG.5x.3{F11.6*10X 

) ) 

116 

FORMAT 

{/17HRA ASC NODE 

DEG.5X*3(F11.6»10X 

) ) 

117 

FORMAT 

(/17HARG PERIGEE 

DEG*5x.3(F11.6.10X 

) ) 

118 

FORMAT 

(/17HMEAN ANOMALY 

DEG*5X*3(F11.6*10X 

) ) 

119 

FORMAT 

(/17HSEMIMAJ AXIS 

ER*7X»3(F11.8*10X 

) ) 

120 

FORMAT 

(/16H J.D.S. UT2W *FU.5) 


121 

FORMAT 

(24X.28HJOHN F. KENNEDY SPACE CENTER) 


122 

FORMAT 

(29X*18HCOMPUTATION BRANCH) 


123 

FORMAT 

(33X.4(A3) ) 



124 

FORMAT 

(26X*8(A3) ) 



125 

FORMAT 

(////) 



127 

FORMAT 

(28X.20HSU8-SATELLITE POINTS) 


129 

FORMAT 

(F3.0) 




FORMAT 

STATEMENTS FOR OUTPUT 


126 

FORMAT 

(56HTIME( J.D.S.) 

LAT(DEG) LON(DEG) 


1 ) ) 




128 

FORMAT 

(/(F12.5*2(7XF7.2) 

1 *6X*F10.1) ) 



C COMPUTE NO. OF TIMES 

POINTS=(T3-Tl)/DELTAT _ 

IF(CNTRL) 500.510.510 
510 WRITE TAPE 7. POINTS* T3. Tl. DELTAT 
500 CONTINU€ _ 

c convert input data and initialize 


HEIGHTIKM 
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RAD=57. 2957795 


CUL=6378«-I65 

SHR=ISHR 

SMM=ISMM 

ST=(SHR + ( (SMM+SSS/60* )/60. ) ) *15, / RAD 

PRINT 105 

PRINT 126 

DO 200 I=l»5 

DO 200 J=1.3 

200 AL( I »J)=EL( I >J) 

T=T1 

LL=POINTS 

L^O 

DO 900 Nsl.LL _ 

C COMPUTE ORBITAL ELEM ENTS AT T 

CALL ELIN(T»T1»T2.T3»AL«SAT>SET>SIT*C0T»S0T) 

AMT=AMINF(T>T1>T2>T3>EL (6»1) . EL ( 6 » 2 ) . EL ( 6 . 3 ) .P 1 ♦ P2 .P3 *D1 . D2 » D3 ) 

C TEST MEAN ANOMALY TO SEE IF ZERO 

IF(AMT) 5»950»5 

5 AMT=AMT/RAD 

C COMPUTE ECCENTRIC ANOMALY AT T 

E T=EKEPF(AMTtSET ) 

SIT = SIT/RAD 

COT=COT/RAD 

SOT=SOT/ RAD 

C COMP UTE INERTIAL COORDINATES AT T 

CALL PV0E(SIT>SAT »SET»S0T»C0T»ET»X»Y»2) 

DELTR=(T-RT)*86400. 
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c 


COMPUTE LAT.» LON. AND H 


CALL GEODT(X»Y»Z»ST» DELTR»VE»ES»ELQN»RHO»6LATl»GLAT2»H»RAfGST) 

ELON = ELON*RAD 

GLAT1=GLAT1*RAD 

GLAT2=GLAT2*RAD _ 

RA = RA»RAD 

GST = GST*RAD 

H = H*CUL 

C PRINT AND/OR WRITE OUTPUT 

IF( CNTRL) A20»405»410 

405 WRITE TAPE 7 » T » GLAT2 > EL ON > H 

T=T+DELTAT 

GO TO 900 

410 WRITE TAPE 7 » T t GLAT 2 f ELON * H 

420 PRINT 126>T»GLAT2»ELONfH 

T = T + DELT AT 

L = L + 2 

IF (L-50) 9Q 0»900»30Q 

300 PRINT 105 

PRINT 126 
L=0 

900 CONTINUE 

950 STOP 

END 


* FORTRAN 

C ELIN-SUBROUTINE TO INTERPOLATE FOR ECCEN.* INCL#»LON. 
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C OF ASC. NODEf ARG. OF PERIGEE. SEMI MAJOR AXIS 

C UTILIZES EVERET FUNCTION SUBPROGRAM 

SUBROUTINE ELI N ( T , T 1 . T2 » T3 , AL .SAT .SET *SI T .COT . SOT ) 

^MENSION AL(5.3).ELT(5) 

DO 5 1=1.5 

5 ELT ( I )=E VERET( T.T1.T2.T3.AL( I .1 ) .AU I .2) .AL( I .3 ) ) 

SAT=ELT(1) 

SET = ELT f2) 

SIT=ELT(3) 

COT= ELT(») 

^T = ELT(5) 

RETURN 

END 


* FORTRAN 

C EVERET-FUNCTION SUBPROGRAM TO INTERPOLATE FOR 1 VALUE 

C GIVEN 3 USIN G EVERETTS 2ND CE NTRAL DIFFERENCE 

C FORMULA 

FUNCTION EVERETIT.Tl. T2.T3. XI. X2.X3) 

W=T2-T1 

5=(T-T2)/W 

R=l.-S 

D=X3-2.* X2+X1 

XI=R*X2+R*( R»R-1. )*D/6.+S*X3+S*(S*S-l. )*D/6. 



EVERET«XI 





RETURN 





END 
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•FORTRAN 


AMIN - FUNCTION SUBPROGRAM FOR 

MEAN ANOMALY INTERPOLATION 


FUNCTION AMIN(T»T1»T2>T3>AM1»AM2>AM3»P1>P2»P3»PD0T1 «PD0T2>PD0T3) 

Cl=518400. 

C2=373*248 

IF (T--T2) 10,20>30 

10 GO TO 80 

20 AMT=AM2 

GO TO 90 

30 IF (T-T3) 40*50»60 

40 GO TO 81 _ 

50 AMT=AM3 

GO TO 90 

60 AMT^'O, 

GO TO 90 

30 0ELT=T-T1 

AMI=AM1 _ 

PI=P1 

PDOTI=PDOTl 

*GO TO 82 ^ 

81_DELT = T-T2 

AMI=AM2 _ 

PI=P2 


PDOTI=PDOT2 



82 AMT=AMI+C1/PI*DELT-C2*PD0TI/(PI*PI )*DELT*DELT 

MODPI=AMT/360. 

FLMOD=MODPI 

AMT=AMT-FLMOD*360. 

90 AMIN = AMT 

C AMI N IS ZERO WHEN T IS GREATER THAN T3 

RETURN 
END 


* FORTRAN 

C EKEP-FUNCTION SUBPROGR AM TO SOLVE KEPLERS E 9 . FOR AN ELLIPSE 

FUNCTION EKEP (AM. ECO 

DIMENSION E( 10 ) 

TOL=. 00000002 

_ Z=ECC*SINF(AM)/SQRTF(ECC»ECC+1.-2.*ECC*C0SF(AM) ) 

C0TAM=C0SF(AM)/SINF(AM) 

E(l)=AM+Z-Z**A*C0JAM/6. 

1 = 1 

10 E( I + l )=E( I ) + i AM+ECC*SINF(E(I ) )-E( I ) [/ ( l.-ECC*COSF(E( I ) ) ) 

DE1JAE=ABSF(E( I+D-E! I ) ) 

IF(DELTAE-TOL) 30,30.20 
20 1 = 1+1 
GO TO 10 
30 EKEP= E(I+i) 

RETURN 

END 
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* FORTRAN 

C PVOE-SUBROUTINE TO COMPUTE POSITION VECTOR FROM ORBITAL ELEMENTS 

C 

SUBROUTINE PVOE ( S I T »SAT ♦ SET »SOT »COT »ET ♦ X » Y *2 ) 

SISIT=SINF(SIT) 

COSIT=COSF(SIT) 

SISOT=SINF(SOT) 

COSOT=COSF(SOT) 

SICOT=SINF(COT) 

COCOT=COSF(COT' 

SIET=SINF(ET) 

COET=COSF(ET) _ 

A=COET-SET 

B=SQRTF( l.-SET»SET)*SIET 

X=SAT*( A*(COSOT*COCOT-SISOT*SICOT*COSIT)+B*(-SISOT*COCOT-COSOT*SK 
lOT^COSIT)) _ 

Y=SAT*< A*(COSOT*SICOT+SISOT*COCOT»COSIT ) +B»(COSOT»COCOT»COSIT-SISO 

1T*SIC0TI) 

Z=SAT*( A*SISOT*SISIT+B*COSOT*SISIT ) 

RETURN 

END 


* FORTRAN 

C 

C GEODT- SUBROUTINE TO COMPUTE LAT., LON, AND HEIGHT 

C FROM inertial X.Y.Z 
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SUBROUTINE GEODT ( X * Y »Z »ST * DT f VE » E2 ♦ ELON > RHO , GL AT 1 #GLAT2 * H * R A , GST ) 


PI=3. 14159265 
TWOPI=6. 28318531 

COMPUTE R,A. OF SATELLITE, 0-360 
rA=ARTNF(Y,X ) 

IF(RA)5,10,10 
5 RA=TWOPI+RA 
GO TO 15 

10 RA = RA .... 

15 CONTINUE 

COMPUTE CORRECTION FOR ROTATION OF EARTH, 0-360 

rt=ve*dt 

IRT=RT/TWOPI 

FRT=IRT 

rt=rt-frt*twopi 

compute longitude, +OR-180 
ELON=RA-ST-RT 
IF(ELON)70,40,90 
70 IF(ELON+PI )80,40,40 
80 ELON=ELON+TWOPI 
GO TO 40 

90 IF(ELON-PI )40,40,30 
30 ELON=ELON-TWOPI 
40 CONTINUE 

COMPUTE GEOCENTRIC LAT., +OR-90 
GLAT1=ATANF(Z/SQRTF(X*X+Y*Y) ) 

COMPUTE GEODETIC LAT. 


E4=E2*E2 



E6 = E2^^E4 


£8=E2*E6 

RHO=SQRTF(X*X+Y*Y+Z*Z ) 

RH02=RH0»RH0 
RH03=RH0*RH02 
RH04- = RH0*RH03 

A2= ( 512.*E2+128.*EA+60.*E6+35.*E8 ) / ( 1024.*RHO) 

1+(E6 + E8 )/( 32,*RH02 ) -3 . * ( 4 . *E6 + 3 . ^i-E8 ) / (256.*RH03 ) 

A4 = -( 64.*E4 + 48.*E6 + 3 5.*E8 ) / ( 1 024 . *RHO ) + ( 4 . ^<-£4+2 . *E6 + E 8 ) / ( 16.*RH02) 
l+( 15.*E8 ) / ( 256.*RH03)-E8/ ( 16.*RH04) 

GLAT2 = GLAT1 + A2*SINF ( 2 . *GL AT 1 ) +A4*S I NF ( 4 . *GLAT 1 ) 

C COMPUTE HEIGHT IN C.U.L. 

S2GLAT = SINF(GLAT2)*SINF(GLAT2 ) 

H=RH0*SINF(GLAT1 ) /SINF(GLAT2)-( 1.-E2 ) /SQR TF ( 1 . -E 2*S2GL AT ) 

RETURN 

END 
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